density.est=function(data,y,g,a,A,units,estmethod){
  data=as.name(data); y=y; g=g; a=a; A=A
  units=as.name(units)
  if(estmethod=='quadrat'){
    n=sum(y); mi=y*g; mbar=sum(mi)/n
    lambda.hat=round(mbar/a,6)
    lambda.vhat=lambda.hat/(a*n)
    lambda.bound=round(2*sqrt(lambda.vhat),6)
    lambda.lower=lambda.hat-lambda.bound; lambda.upper=lambda.hat+lambda.bound
    Mhat=ceiling(lambda.hat*A)
    Mvhat=A^2*lambda.vhat
    Mbound=ceiling(2*sqrt(Mvhat))
    Mlower=Mhat-Mbound; Mupper=Mhat+Mbound
  } else if(estmethod=='stockedq'){
    if(length(g)==1){n=g}
    else if(length(g)>1){n=length(g)}
    if(length(y)==1){y=y}
    else if(length(y)>1){y=sum(y)}
    lambda.hat=-(1/a)*log(y/n)
    lambda.vhat=(1/(n*a^2))*(exp(lambda.hat*a)-1)
    lambda.bound=2*sqrt(lambda.vhat)
    lambda.lower=lambda.hat-lambda.bound; lambda.upper=lambda.hat+lambda.bound
    Mhat=ceiling(lambda.hat*A)
    Mvhat=A^2*lambda.vhat
    Mbound=ceiling(2*sqrt(Mvhat))
    Mlower=Mhat-Mbound; Mupper=Mhat+Mbound
  }
  cat("","\n","Population estimation results of",estmethod,"method",'\n',"Data:",data,"\n","Per unit area =",a,units,
      "Total area =",A,units,"\n","Density estimate =",lambda.hat,'\n',"Vhat density =",lambda.vhat,"\n",
      "Bound (density) =",lambda.bound,"\n","Lower Bound (density) =",lambda.lower,"Upper Bound (density) =",
      lambda.upper,"\n","Total number of elements =",Mhat,'\n',"Vhat total =",Mvhat,"\n","Bound (total) =",Mbound,
      "\n","Lower Bound (total) =",Mlower,"Upper Bound (total) =",Mupper,"\n","")
  results=list(estmethod=estmethod,data=data,units=units,lambda.hat=lambda.hat,lambda.vhat=lambda.vhat,
               ambda.bound=lambda.bound,lambda.lower=lambda.lower,lambda.upper=lambda.upper,Mhat=Mhat,
               Mvhat=Mvhat,Mbound=Mbound,Mlower=Mlower,Mupper=Mupper)
}
# to use the function with its call:
# density.est(data,y,g,a,A,units,estmethod)

# data: name of dataset, in quotes
# y: vector of frequencies or count of successes (single value)
# g: vector of quadrats or total number of quadrats (single value)
# a: per unit area, single value
# A: total area, single value
# units: units of measurement, in quotes
# estmethod: c('quadrat','stockedq')